/*
	Линейная интерполяция
*/

#include <math.h>
#include <assert.h>
#include <stdlib.h>
#include "linpol.h"

/* Создание объекта полинома */
linpol_t *linpol_new(const int linpol_degr)
{
	assert(linpol_degr > 0);

	/* Выделение памяти под полином */
	linpol_t *linpol_n = calloc(1, LINPOL_T_SZ);
	if(!linpol_n) return 0;

	/* Выделение памяти под данные полинома */
	void *mem_n = 0;

	/* Выделение памяти под значения X */
	mem_n = linpol_alloc_mem(linpol_degr + 1);
	if(!mem_n) return 0;

	linpol_n->vals_x = mem_n;

	/* Выделение памяти под значения Y */
	mem_n = linpol_alloc_mem(linpol_degr + 1);
	if(!mem_n) return 0;

	linpol_n->vals_y = mem_n;

	/* Установка параметров полинома */
	linpol_n->linpol_degr = linpol_degr;

	return linpol_n;
}

/* Удаление объекта полинома */
void linpol_destroy(linpol_t *linpol_o)
{
	assert(linpol_o);

	free(linpol_o->vals_x);
	free(linpol_o->vals_y);
	free(linpol_o);
}

/* Генерирование исходных данных */
void linpol_gen_data(const linpol_t *linpol_o, double arg_x_start, double arg_x_end, double (*gen_func) (const double arg_x))
{
	assert(linpol_o);
	assert(gen_func);

	/* Количество исходных точек */
	int pnts_c = (linpol_o->linpol_degr);

	/* Интервал значений */
	double arg_x_step = (arg_x_end - arg_x_start) / (double) pnts_c;

	/* Генерирование данных */
	for(int i = 0; i <= pnts_c; i++)
	{
		*((linpol_o->vals_x) + i) = arg_x_start;
		*((linpol_o->vals_y) + i) = gen_func(arg_x_start);

		arg_x_start += arg_x_step;
	}
}

/* Вычисление значения полинома */
double linpol_get_val(const linpol_t *linpol_o, const double arg_x)
{
	assert(linpol_o);

	/* Количество рабочих точек */
	int pnts_c = (linpol_o->linpol_degr) + 1;

	/* Выделение памяти под треугольную матрицу значений */
	double **linpol_vals = calloc(pnts_c, sizeof(double *));
	assert(linpol_vals);

	*linpol_vals = linpol_o->vals_y;

	for(int i = 1; i < pnts_c; i++)
	{
		*(linpol_vals + i) = linpol_alloc_mem(pnts_c - i);
		assert(*(linpol_vals + i));
	}

	/* Вычисление значения полинома */
	double linpol_res_val = *(linpol_o->vals_y);
	double linpol_res_tmp = 1;

	for(int i = 1; i < pnts_c; i++)
	{
		for(int j = 0; j < pnts_c - i; j++)
		{
			/* Разность Y */
			double val_y = linpol_vals[i - 1][j + 1] - linpol_vals[i - 1][j];

			/* Разность X */
			double val_x = linpol_o->vals_x[i + j] - linpol_o->vals_x[j];

			/* Результат */
			linpol_vals[i][j] = val_y / val_x;
		}

		linpol_res_tmp *= arg_x - linpol_o->vals_x[i - 1];
		linpol_res_val += **(linpol_vals + i) * linpol_res_tmp;
	}

	/* Освобождение памяти */
	for(int i = 1; i < pnts_c; i++)
	{
		free(*(linpol_vals + i));
	}

	free(linpol_vals);

	/* Возврат результата */
	return linpol_res_val;
}

/* Вычисление теоретической погрешности */
double linpol_get_tfault(linpol_t *linpol_o, const double arg_x, double (*dgen_func) (const double arg_x))
{
	assert(linpol_o);

	/* Количество исходных точек */
	int linpol_degr = linpol_o->linpol_degr + 1;
	double deriv_max = dgen_func(*(linpol_o->vals_x));

	/* Поиск максимума производной */
	for(int i = 0; i < linpol_degr; i++)
	{
		double deriv_t = dgen_func(*(linpol_o->vals_x + i));
		deriv_max = (deriv_t > deriv_max) ? deriv_t : deriv_max;
	}

	/* Подсчет факториала */
	int linpol_fact = 1;

	for(int i = 1; i <= (linpol_degr + 1); i++)
	{
		linpol_fact *= i;
	}

	/* Вычисление погрешности */
	double linpol_fault = 1;

	for(int i = 0; i < linpol_degr; i++)
	{
		linpol_fault *= (arg_x - *(linpol_o->vals_x + i));
	}

	return deriv_max * linpol_fault / linpol_fact;
}